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Abstract 

The emergence and vanishing of superdiffusion in quasi-two-dimensional Yukawa systems are 
investigated by molecular dynamics simulations. Using both the asymptotic behaviour of the 
mean-squared displacement of the particles and the long-time tail of the velocity autocorrelation 
function as indicators for superdiffusion, we confirm the existence of a transition from normal 
diffusion to superdiffusion in systems changing from a three-dimensional to a two-dimensional 
character. A connection between superdiffusion and dimensionality is established by the behaviour 
of the projected pair distribution function. 

PACS numbers: 52.27.Gr, 52.27.Lw, 82.70.Dd, 66.10.cg 
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I. INTRODUCTION 



Diffusion is a fundamental transport mechanism which plays a dominant role in many 
physical, chemical and biological systems. It is not only of academic but also practical 
interest to study diffusion in two-dimensional systems since many real-world systems can be 
described as being two-dimensional or quasi-two-dimensional, including surfaces or layers 
of small width, e.g. quantum wells. Two-dimensional diffusion has long been known to 
exhibit anomalous behaviour for a number of interactions and systems. One of the possible 
anomalies is the so-called super diffusion which describes diffusion proceeding faster than 
normal diffusion in the sense that a particle achieves a greater distance from its starting point 
than expected from Pick's law. For three-dimensional simple systems in thermodynamical 
equilibrium, to the best of our knowledge, no super diffusion has been observed to date. On 
the other hand, under nonequilibrium conditions, anomalous transport is well-known, e.g. 

n n 

for chaotic systems [11], turbulent flows [2f], or plasmas in turbulent magnetic fields [3]. These 
will not be considered here. For systems which are neither three-dimensional nor strictly 
two-dimensional, we expect a gradual transition from super diffusive behaviour to Fickian 
diffusion. 

In two-dimensional systems of hard disks, a slow oc decay of the long-time tail of 
;he velocity autocorrelation function (VACF) was first observed by Alder and Wainwright 
^. Such a decay results in a divergent Green-Kubo integral. Superdiffusion was found to 
take place in quasi- 

systems [oj] and systems under laser-induced shear |lO|. No experimental evidence for su- 
perdiffusion was found by Nunomura et al. who examined an underdamped liquid complex 
plasma 



The intent of this study is to examine the superdiffusion in systems which are quasi- 
two-dimensional, i.e. the extension in one spatial dimension of the system is much smaller 
than in the other two. This is often the situation in experimental setups, for example in 
dusty plasma experiments where the dust grains are levitated by an electrostatic force which 
is counteracted by gravity. It is clear that under such circumstances the particles are not 
rigorously restricted to a two-dimensional plane but form a quasi-two-dimensional system. 
For example, Ref. 6|] reports on experiments in which the dusty plasma under consideration 
consisted of two to three layers. 
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Here we consider a macroscopic system of charged particles interacting via a screened 
Coulomb (Yukawa) potential. It is of relevance to dusty plasmas [12], colloidal suspensions, 
electrolytes and other systems. It is also an important theoretical tool since it allows to 
tune the inter-particle interaction (by varying the screening length of the Yukawa potential) 
from being very long-ranged to almost contact interaction. 

We report on molecular dynamics studies performed for quasi-two-dimensional systems. 
Our focus lies on the influence of the degree of quasi-two-dimensionality on the vanishing 
and emergence of superdiffusion. 



II. MODEL AND SIMULATION TECHNIQUE 

We study the system using equilibrium molecular dynamics simulations (e.g. [13|]). The 
interaction of the particles is given by the Yukawa pair potential 

0(r) = (1) 

Here, Q is the particles' charge and Ad is the screening length. 

A two-dimensional Yukawa system is characterized by two parameters, the Coulomb coupling 
parameter V = {Q"^ / Atteo) x (l/a^skBT) and the screening parameter k = a^s/^D- T is the 
temperature and a^s = (nvr)"^/^ is the Wigner-Seitz radius for two-dimensional systems. 
While our simulation box is periodic in the x and y directions, it is unbound in the z 
direction. The particles' perpendicular movement in z direction is restricted by one of two 
confinement potentials 

r)2 /p-{z+w)/XD f,-{-z+w)/\D\ 
V'-iz) = ii- "—-^ + (3) 



47r£o \ z + W -z + W ^ 

The confinement ([2]) is a simple harmonic trap where / denotes the trap amplitude jl^ . 
The second model, Eq. ([3]), uses a box-shaped confinement with "soft" walls. It can be 
thought of as consisting of two rigid walls of immovable particles with the same interaction 
as the actual particles, i.e. a Yukawa interaction, see Fig. [1] for an illustration. Here, W 
is defined as W = ^ + a^s with w being the "width" of the box, i.e. twice the maximum 
displacement in +z or —z direction. We use the width w as a third parameter for the 
soft-box confinement. 
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FIG. 1: Soft-box potential (piz) and corresponding particle density n(z). The trap is impenetrable 
at ±W and the particles can typically explore the trap width w (here, w = 2a^s). 

The choice between two confinements allows us to (a) reproduce typical experimental 
situations or (b) study the effect of the dimensionality on superdiffusion while retaining 
comparable plasma conditions at all times. We achieve this by the following procedure: In 
the harmonic confinement, as the amplitude of the trap is decreased, the particle number 
is left unchanged and the particles can explore a wider vertical space. This results in an 
increased mean inter-particle distance, i.e. the (three-dimensional) density is decreased. 
Contrarily, in the soft-box confinement, we chose point to fix for the 3D 

density and scale the particle number as necessary to maintain that density for all values of w. 




We performed simulations for a fixed coupling parameter T = 200. Prior to measure- 
ment, the particles' velocities are rescaled at each timestep to the desired temperature until 
a Maxwell distribution is well-established. During the measurement, the velocities are not 
rescaled. Our simulations are carried out for k = 2.0 and k = 3.0. At these parameters, a 
2D Yukawa system is well in the liquid phase, with the melting points being F ^ 415 and 
F ^ 1210 for K = 2.0 and 3.0 respectively ll5|. 

The particle number in our simulations is = 6000 for the harmonic confinement and 

= 4000 . . . 16000 for the soft-box confinement with A^ > 6000 for w > 1.5a^,. 
In the following, time is given in units of the inverse plasma frequency ujp = 
{Q"^ 127X8 Qma^^Y^'^ with m being the mass of the particles and lenghts are measured in 
units of a ,„.9. We solve the equation of motion for each particle using the velocity Verlet 
algorithm 16|. 

The structure of the systems is characterized by measuring the density distribution n{z) 
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perpendicular to the confined direction and by the projected pair distribution function g*{r). 
To analyze diffusion properties we calculate the time-dependence of the mean-squared dis- 
placement (MSD) 

Ur{t) = im - f{to)\') (4) 

where (.) denotes an ensemble- average and r = {x y) is the position vector in the plane; the 
z-component is bound from above due to the confining force and does not need to be taken 
into account. 

The motion can be classified according to the time-dependence Ur{t) ~ f^. Normal Fickian 
diffusion is characterized by a linear time dependence, 7 = 1. If7>lor7<l, motion 
is super- or subdiffusive, respectively. Ballistic, i.e. undisturbed, motion is trivially marked 
by 7 = 2. 

We have carried out calculations of the MSD for different trap amplitudes and box widths 
and determined the slope of the curve on a double logarithmic plot between t = lOOa;"^ and 
t = 300ci;~^ which yields the diffusion exponent 7. 

A more direct insight into a particle's movement can be obtained from the decay of the 
velocity autocorrelation function (VACF) 

Z{t) = {v{t) ■ vito)) (5) 

where v = {v^ Vy). Z{t) is a measure of the memory of the system. For uncorrelated binary 

collisions, Z{t) is expected to decay exponentially. If Z(t) decays algebraically, Zit) ~ 

a needs to be larger than 1 for diffusion to be Fickian, because otherwise no valid diffusion 



coefficient can be calculated from the Green-Kubo formula. For 3D, the decay has been 
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to be algebraic with a = 1.5 for a number of pair potentials in experiments llSl, Il9| . 
simulations and by theoretical models [2l|, |22[ |. 



For 2D Yukawa systems exhibiting super diffusion, Liu and Goree have found an algebraic 
decay with a ^ 1 [23|. This is an indication of super diffusive behaviour. 

In our simulations, we calculated the velocity autocorrelation for different trap ampli- 
tudes. To obtain accurate statistics, we performed between 10 and 50 runs of 700.000 
timesteps for each trap amplitude and again determined the slope of the curve in a double 
logarithmic plot, this time between t = lOOcu"^ and t = 2500;"^ Special attention must be 
paid to the error estimation. We used the Jackknife method to evaluate our data, because 



standard error estimates may not be sufficient in this case 
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FIG. 2: The MSD over time for different trap amplitudes in tlie harmonic confinement and k = 2.0. 
The sohd hnes have slope 2.0 (ballistic regime) and 1.0 (normal diffusion) respectively. The slope 
of the MSD in this log-log plot is the diffusion exponent, see Fig. [3l 

III. RESULTS 

A. Harmonic confinement 

We begin by noting that particles interacting via a Yukawa potential and confined by 
a harmonic trap support the formation of layers [3]. The number of layers formed 
depends on the temperature, the trap frequency and the screening length. For our choice of 
parameters, we found that for n = 2.0 a second layer builds up at / = 0.12 while for k = 3.0 
the system consists of a single layer for the whole range of / examined (cf. the insets above 
Figs. [3]and|lD. 

Typical results for the time-dependence of Ur{t) are shown in Fig. [2l In this double- 
logarithmic plot, the slope of the curves corresponds to the exponent of the algebraic be- 
haviour. For if: < lOujp^, Ur{t) grows quadratically with time, which corresponds to a ballistic 
motion of each individual particle. After a narrow transition region of about 10uj~^, the par- 
ticles' movement is dominated by diffusive processes. The slope in this region allows us to 
classify motion as superdiffusive, diffusive or subdiffusive, i.e. it is the diffusion exponent 7. 
In Fig. [21 Ur{t) is depicted for different trap amplitudes /. For strong confinements, Ur{t) 
deviates strongly from a purely diffusive behaviour of 7 = 1. For increasingly more relaxed 
confinements, the slope of Ur{t) tends more and more towards unity, that is the migration 
becomes less superdiffusive. 




6 




+ + + + + + 




trap amplitude / 



FIG. 3: The diffusion exponent for different trap amplitudes and k = 2.0 in tlie harmonic con- 
finement. The straight line is a guide for the eyes. Also shown is the value of the (projected) 
pair correlation function at the distance r = a^js- The top graphs show the density profile in the 
confined direction from z = — 4a^s to z = 4:aws at the trap amplitude indicated by the arrows. 
The n{z) distributions are normalized here to unit amplitude. 

The dependence of the diffusion exponent 7 on the trap amplitude / is shown in Figs. [3] 
andU 

It is clear that the degree of superdiffusivity, as measured by the diffusion exponent 
7, decreases with an increasing width of the system. For k, = 2.0, Fig. [3l the diffusion 
exponent is in the vicinity of 7 = 1.16 for strong confinements which coincides with the 
value we obtained for strictly 2D systems. When the width of the particle distribution along 
the z-axis exceeds 2a^s (which is near the mean interparticle separation in the strictly 2D 
case), the diffusion exponent begins to deviate from that value and continues to fall for lower 
trap amplitudes /. The same behaviour can be seen in Fig. IHfor k = 3.0. Here, 7 saturates 
at ~ 1.20 for strongly confined systems and again falls when the system width exceeds 2a^s. 
As noted before, the 3D density of the system differs for different trap amplitudes. To 
exclude the possibility that the vanishing of superdiffusion is only an effect of the density, 
we also simulated strictly 2D systems in which we changed the particle density until the first 
peak in the pair distribution function coincided with that of our quasi-2D systems. Our data 
(not shown) clearly indicates that a decreased density does not lower the diffusion exponent 
in our parameter range. In fact, for k = 2.0 the superdiffusion is stronger for systems of 
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FIG. 4: The diffusion exponent for different trap amplitudes and k = 3.0 in the harmonic confine- 
ment. Also shown is the value of the (projected) pair correlation function at the distance r = aws- 
The top graphs show the density profile from z = — 4a^s to z = 4a^s in the confined direction at 
the trap amplitude indicated by the arrows. The n{z) distributions are normalized here to unit 
amplitude. 

lower density. 

To support the idea that the vanishing of super diffusion is connected with the fact that 
particles can pass each other in z-direction, we analyze the projected pair distribution func- 
tion g*{r). Its value at r = a^s is also shown in Figs. [3]and|H For nearly 2D systems (high 
/), g*{aujs) is practically zero. The reason for this is that for strongly correlated liquids, the 
particles' kinetic energy is not sufficient for two particles to come together as close as one 
Wigner-Seitz radius. Instead, they are trapped in local potential minima from which they 
can escape only after some time. Figs. [3] and H] show that g*{au]s) is non-zero for lower trap 
amplitudes / < 0.1. This is a result of the projection of all particles onto a two-dimensional 
plane. Individual particles are still separated by more than a^s except now the finite width 
of the system allows the particles to go "over and under" each other in z-direction and the 
particles' projections can be close. 

By inspection of Figs. [3]and|U we see that the behaviour of the dynamic quantity 7(/) 
is closely mirrored by the static quantity g*{aws)\f- This underlines that superdiffusion is 
connected with the dimensionality of the system. 



Fig. 5(a) depicts the asymptotic behaviour of Z{t) for different trap amplitudes / in a 
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FIG. 5: [(a)] Double-logarithmic plot of the long-time tail of the VACF Z{t) (Z(0) = 1) from 
/ = 0.01 (top) to / = 0.28 (bottom). 

|(b)| The exponent of the algrebraic decay of the velocity auto-correlation function for k = 2.0 and 
different trap amplitudes in the harmonic confinement. Error bars denote standard errors from 
the Jackknife estimator. The top graphs show the density profile in the confined direction from 
z = — 4a^s to z = 4a^s at the trap amplitude indicated by the arrows. The n{z) distributions are 
normalized here to unit amplitude. 



double logarithmic plot. The results for the decay of the VACF are shown in Figs. 5(b) and 
EioT K = 2.0 and k = 3.0. 



The data shown in Fig. 5(a) validate our approach of modelling the asymptotic behaviour 



of Z{t) as an algebraic decay since Z{t) closely follows a straight line in the log-log plot. 



The peak seen in Fig. 5(a) at long times for high / is due to our finite simulation box and 



is caused by sound waves travelling through and re-entering the system due to the periodic 
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FIG. 6: The exponent of the algrebraic decay of the velocity auto-correlation function for k = 3.0 
and different trap amplitudes in the harmonic confinement. Error bars denote standard errors from 
the Jackknife estimator. The top graphs show the density profile in the confined direction from 
z = —iaws to z = 4a^s at the trap amplitude indicated by the arrows. The n{z) distributions are 
normalized here to unit amplitude. 



boundary conditions. Measurement is limited to times smaller than the time of the sound 



wave traversal. The dependence of a on the trap amplitude /, Fig. 5(b) and El indicates 
vanishing of super diffusion for broader systems. Starting close to the value a = 1.0 for 
narrow systems, we see an increase in a for broader systems. This corresponds to a faster 
decay which is indicative of a loss of the particles memory, i.e. at subsequent times, a particle 
is less likely to travel in its original direction. 



B. Soft-box confinement 



We now turn our attention to the case of the soft-box confinement. Again, let us note 
that here too, the system forms layers when given enough space in the confined direction. 
We find that the number of layers formed is higher than in the harmonic confinement which 
we attribute to the constancy of the particle density. Recall that in this case we change the 
number of particles to maintain a constant 3D density. 

The dependence of the diffusion on the width of the system. Figs. [7] and [H is more 
involved than for the harmonic confinement. Again, the general trend is for superdiffusion 
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FIG. 7: (color online) The diffusion exponent for different box widths and k = 3.0 in the soft-box 
confinement. Also shown is the value of the (projected) pair correlation function at the distance 
1^ = o,ws- The top graphs show the density profile in the confined direction from z = — Sa^u^ to 
z = 3a^s at the box width indicated by the arrows. The n{z) distributions are normalized here to 
unit amplitude. The background color separates regions of one, two and three layers. 



to vanish for increasingly broader systems. But here, the vanishing happens in stages: After 
a first drop, the diffusion exponent reaches a plateau from which it drops to a second plateau. 
This behaviour is connected to the formation of layers in the system as indicated by the 
different background colors in Figs. [7] and [8] (cf . top graphs in these Figs). As another layer 
is formed in the system, the diffusion exponent experiences a drop. The non-monotonicity 
of the curve in Figs. [7| and [8] is due to statistical error. 

We see that two to three layers are already sufficient to reduce the superdiffusive behaviour 
substantially. In addition, we again notice that the value of the projected pair distribution 
function g*{au,s) is directly correlated with the diffusion exponent. 



IV. SUMMARY AND DISCUSSION 



A study of superdiffusion in quasi-two-dimensional Yukawa liquids was performed by equi- 
librium molecular dynamics simulations. The two indicators for superdiffusion employed, 
the MSD and the VACF, both show sensivity to the dimensionality of the system. For 
increasingly broader systems superdiffusion gradually vanishes. The transition from su- 



11 



1.20 



§ 1.15 

o 
ft 

X 

« 1.10 

o 

i 1.05 

^3 



1.00 




1.5 2 2.5 3 3.5 4 

box width w/ams 



FIG. 8: (color online) The diffusion exponent for different box widths and k = 2.0 in the soft-box 
confinement. The top graphs show the density profile in the confined direction from z = —Sa^s to 
at the box width indicated by the arrows. The n{z) distributions are normalized here to 
unit amplitude. The background color separates regions of one, two and three layers. 



perdiffusion to normal diffusion was tested for two representative values of n and found to 
be qualitatively comparable. This leads us to the conclusion that the transition is universal 
for Yukawa systems in the fluid phase. 

To ensure that the choice of confinement does not interfere with the change in dimen- 
sionality, we have used two different schemes to confine the system. The general trend of 
the vanishing of super diffusion appears to be independent of the type of confinement. The 
finer details of how the vanishing happens depend on the choice of confinement and here 
especially on the formation of layers. 

The strength of superdiffusion at zero width and the number of layers in the systems 
depend on the inter-particle potential and the system parameters. At a fixed Coulomb 
coupling parameter T = 200, superdiffusion is stronger for k, = 3.0 and weaker for k, = 2.0. 

By inspection of the projected pair distribution function we have established a close 
connection between superdiffusion and the dimensionality of the system. To this end, we 
display in Fig. [9] the dependence of the diffusion exponent 7 on the value of the projected 
pair distribution g*{aws) for k, = 2.0 and k, = 3.0 in the harmonic confinement, which is close 
to linear. For the soft-box confinement the general behaviour is similiar, although details 
are more complex due to the occurence of plateaus in the curves, cf. Fig. [71 
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FIG. 9: The diffusion exponent 7 over the value of the projected pair distribution g*{aws) for 
K = 2.0 and k = 3.0 in the harmonic confinement. 

It remains an interesting question what types of other pair potentials also support su- 
perdiffusion and how its strength depends on the interaction range. Finally, it will also be 
of high interest for future analysis to see how quantum effects influence super diffusion. This 
could be done, e.g., by use of effective quantum pair potentials 



27 



28| 



Acknowledgment s 



We acknowledge stimulating discussions with J.W. Dufty. This work has been supported 
by the Deutsche Forschungsgemeinschaft via SFB-TR24 grant A5 and by the Hungarian 
Scientific Research Fund, through grants OTKA-T-48389, OTKA-IN-69892 and PD-049991. 



[1] G. Zaslavsky, Phys. Rep. 371, 461 (2002). 

[2] T. Hauff, F. Jenko, and S. Eule, Phys. Plas. 14, 102316 (2007). 

[3] P. Pommois, G. Zimbardo, and P. Veltri, Phys. Plas. 14, 012311 (2007), and references therein. 

[4] B. Alder and T. Wainwright, Phys. Rev. A 1, 18 (1970). 

[5] S. Ratynskaia, K. Rypdal, C. Knapek, S. Khrapak, A. V. Milovanov, A. Ivlev, J. J. Rasmussen, 

and G. E. Morfill, Phys. Rev. Lett. 96, 105010 (2006). 

[6] R. A. Quinn and J. Goree, Phys. Rev. Lett. 88, 195001 (2002). 



13 



[7] Y.-J. Lai and L. I, Phys. Rev. Lett. 89, 155002 (2002). 
[8] W.-T. Juan and L. 1, Phys. Rev. Lett. 80, 3073 (1998). 
[9] B. Liu and J. Goree, Phys. Rev. Lett. 100, 055003 (2008). 
[10] W.-T. Juan, M.-H. Chen, and L. I, Phys. Rev. E 64, 016402 (2001). 

[11] S. Nunomura, D. Samsonov, S. Zhdanov, and G. Morfin, Phys. Rev. Lett. 96, 015003 (2006). 
[12] M. Bonitz, D. Block, O. Arp, V. Golubnychiy, H. Baumgartner, P. Ludwig, A. Piel, and 

A. Fihnov, Phys. Rev. Lett. 96, 75001 (2006). 
[13] V. Golubnychiy, M. Bonitz, D. Kremp, and M. Schlanges, Phys. Rev. E 64, 16409 (2001). 
[14] Z. Donko, P. Hartmann, and G. Kalman, Phys. Rev. E 69, 65401 (2004). 
[15] P. Hartmann, G. Kalman, Z. Donko, and K. Kutasi, Phys. Rev. E 72, 26409 (2005). 
[16] M. Allen and D. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1987). 
[17] G. Paul and P. Pusey, J. Phys. A: Math. Gen. 14, 3301 (1981). 

[18] C. Morkel, C. Gronemeyer, W. Glaser, and J. Bosse, Phys. Rev. Lett. 58, 1873 (1987). 
[19] J. X. Zhu, D. J. Durian, J. Miiller, D. A. Weitz, and D. J. Pine, Phys. Rev. Lett. 68, 2559 
(1992). 

[20] A. McDonough, S. P. Russo, and I. K. Snook, Phys. Rev. E 63, 026109 (2001). 
[21] Y. Pomeau and P. Resibois, Phys. Rep. 19, 63 (1975). 

[22] R. Zwanzig and M. Bixon, Phys. Rev. A 2, 2005 (1970). 
[23] B. Liu and J. Goree, Phys. Rev. E 75, 016405 (2007). 

[24] Given M samples {yo{t) , yi{t) , . . . yM{t)} of the VACF, we perform M least-square fits of the 
average of M — 1 of these samples to the target function at + c, leaving one different sample 
out from the average each time. The mean value a is computed from the M fit values and the 
error is estimated as the square root of 0-2 = - 1) - af/M where a is the result 

obtained by fitting over the average of all M samples. 

[25] J. Shao and D. Tu, The Jackknife and Bootstrap (Springer, New York, 1995). 

[26] H. Totsuji, T. Kishimoto, and C. Totsuji, Phys. Rev. Lett. 78, 3113 (1997). 

[27] A. Filinov, M. Bonitz, and W. Ebeling, J. Phys. A: Math. Gen. 36 (2003). 

[28] A. Filinov, V. Golubnychiy, M. Bonitz, W. Ebeling, and J. Dufty, Phys. Rev. E 70, 46411 
(2004). 



14 



